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ABSTRACT 

Debris disks should not be completely gas-free, since there is second generation gas from out- 
gassing of planetesimals and dust grains via sublimation, photodesorption, or collisions, generating 
a system of dust-to-gas ratio close to unity, where hydrodynamics cannot be ignored. A clumping 
instability exists in this configuration, that has been hitherto explored only in one-dimensional, in- 
compressible models. We performed 2D numerical compressible models of a disk with comparable 
amounts of gas and dust to study the growth and development of this instability. Our model solves 
the momentum equation for the gas and dust, together with energy and continuity equations. We 
uncover that the backreaction of the drag force from the gas onto the dust shepherds rings, similar 
to those observed in debris disks and usually attributed to the presence of h5^othetical undetected 
planets. We also uncover that the eccentricity of these rings, usually presented as convincing evi- 
dence for the presence of a planet, can actually be simply explained by a standing wave propagating 
along the ring. The rings support a spectrum of oscillations, with one particular mode representing 
epicyclic motion. The apparent eccentricity matches the eccentricity in observed systems. This sug- 
gests that the planet possibility, though thrilling, is not necessarily required to explain these systems. 



1. INTRODUCTION 

Disks around young stars appear to pass through an 
evolutionary phase when the disk is optically-thin and 
the dust to gas ratio is of order unity give or take an or- 
der of magnitude. It can be hard to precisely measure 
the total masses of the dust and gas in such disks, but 
the nearby stars /? Pictoris (Lagrange et al. 1998; Olofs- 
son et al. 2001; Brandeker et al. 2004; Roberge et al. 2006; 
Troutman et al. 2011), HD32297 (Redfield 2007), 49 Ceti 
(Zuckerman et al. 1995) and HD 21997 (Moor et al. 2011) 
all host disks of dust resembling ordinary debris disks 
and also have stable circuutnstellar gas detected in molec- 
ular CO, Na I or other metal lines; the inferred mass of 
gas ranges from Lunar masses to a few Earth masses. 
The gas in these disks is thought to be produced by 
planetesimals or dust grains themselves, via sublima- 
tion or photodesorption (Gregorieva et al. 2007) or colli- 
sions (Czechowski & Mann 2007), processes that should 
occur in every debris disk at some level. 

This kind of disk may form rings, clumps, or spiral 
structures via a recently proposed instability (Klahr & 
Lin 2005; Besla & Wu 2007) that is the subject of this pa- 
per. Under a wide range of conditions, gas drag causes 
dust in a gas disk to concentrate at pressure maxima in 
the disk (Takeuchi & Art5anowicz 2001). But when the 
disk is optically-thin to starlight, the gas is most likely 
primarily heated by the dust, via photoelectric heating. 
In this circumstance, a concentration of dust that heats 
the gas creates a local pressure maximum that in turn 
can cause the dust to concentiate more. 
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The result of this instability could be that the dust and 
gas clump into rings or spiral patterns or other struc- 
tures that could be detected via coronographic imag- 
ing or other methods. Indeed, images of debris disks 
and transitional disks show a range of asymmetries 
and other structures that beg for explanation. Klahr 
& Lin (2005) raised the possibility that the instability 
they hj^othesized could explain some of the observed 
structures. Alternative explanations for these structures 
sometimes rely on planetary perturbers-a tantalizing 
possibility. But we are interested in investigating any 
possible explanation for these disk structures that does 
not require a hidden planetary companion. 

More investigation of this instability is needed. Pre- 
vious investigations of the instability neglected a cru- 
cial aspect of the disk d5mamics: the momentum equa- 
tion for the dust and gas. In other words, they effec- 
tively ignored the inertia of the dust and gas, a property 
that could damp the growth of the instability, or other- 
wise fundamentally alter it. Moreover, prior investiga- 
tions only considered one-dimensional models, which 
can only investigate azimuthally-symmetiic, ring-like 
patterns. This limitation also left open the possibility 
that in higher dimensions, the power in the instability 
might collect in higher azimuthal wavenumbers, gen- 
erating only small clumps that would be imobservable. 
Such behavior would make this instability uninteresting 
as an explanation for the structures seen in circumstellar 
disks. 

To deepen our understanding of this instability, we 
performed a 2-D numerical and analytic investigation of 
a disk containing comparable masses of gas and dust. 
Our model solves the momentum equation for the gas 
and dust, together with energy and continuity equa- 
tions, an improvement over prior models. We report on 
the parameter space where this stability grows (dust-to- 
gas ratio, thermal and coupling times). Finally, we dis- 
cuss the possibility that this instability could result in 
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detectable structures in transitional disks or debris disks 
with gas. 

2. THE PHYSICAL REGIME 

Debris disks with gas represent a regime of nebular 
astrophysics that has only recently been quantified. The 
archetype of the class, the best studied object so far, is 
the disk around j6 Pictoris. Zagorovsky et al. (2010) have 
assembled a model for the properties of this disk that we 
will use as a reference point for our study. 

It is important to note that the total mass of gas that 
these debris disks have is poorly known, even for the 
well-studied disk around /5 Pictoris. Debris disk gas 
has mainly been observed in emission lines from metal 
ions and CO, but the bulk of the gas is generally as- 
sumed to be hydrogen, a component that is difficult to 
measure (See e.g. Thi et al. 2001). Estimates of the hy- 
drogen abundance in the /5 Pictoris disk relative to the 
solar value range from 10"'' to 1; this range translates 
into a range of total gas mass from about 8 x lO^^M^ to 
O.8M0. 

The dust mass is better constrained. Still, care must 
be taken to specify the particle size range of interest, 
since in typical grain size distributions, the larger bodies 
contain most of the mass, yet the smaller grains are the 
ones that we can detect. Zagorovsky et al. (2010) quote 
a dust mass of O.27M0 for particles smaller than 1 cm, 
assuming spherical grains with a number distribution 
dn/ds cx s"^-^, where s is the grain radius. Let us call the 
dust-to-gas ratio e. Given these numbers, e for /5 Pictoris 
would lie in the range of 0.3 to 300. 

These numbers roughly span the range of parame- 
ters for other debris disks with gas, with 49 Ceti a no- 
tably gas-rich exception. (Moor et al. 2011) estimate that 
the HD 21997 disk has a dust mass of ~ O.IM^ and a 
gas mass of O.35M0, corresponding to e ~ 0.3. (Red- 
field 2007) estimates an upper limit to the gas mass for 
HD 32297 of - O.3M0; (Maness 2008) find that the dust 
mass in this system is in the range of 0.02-1 M^, yield- 
ing e in the range ~0.05-3. The 49 Ceti debris disk has 
about the same mass of dust as HD21997, but substan- 
tially more gas: 13M0 (Zuckerman et al. 1995), corre- 
sponding to £ w 0.008. 

Another important quantity for the instability we in- 
vestigated is the thermal time scale. In the Zagorovsky 
et al. (2010) model (the "j6 Pictoris" model, not the "fidu- 
cial" model), the dust is concentrated in a ring about 
100-140 AU from the star. At the peak in the dust den- 
sity, the midplane gas density is about 10 cm~^, the 
dust temperature is roughly 100 K, the gas tempera- 
ture is roughly 70 K. The gas is primarily heated by 
photoelectric emission from dust grains, and primarily 
cooled through the C II 157.7 j.im line emission and the 

total heating/ cooling power is roughly 2 x 10"^^ erg s~^ 
cm"^. Since the specific heat of molecular hydrogen at 
70 K is roughly 1.3 x 10^ erg g"^ K~^, the thermal time 
scale in this model is about 0.5 (n / lOcm"'') years, where 
n is the midplane H2 number density. Given the range 
of possible hydrogen abimdances in the disk, the range 
of time scales of interest corresponds to about 10"^ to 0.1 
orbital periods. 
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Fig. 1. — Two of the five solutions of Eq. (32), corresponding to 
damped oscillations through most of the parameter space. In a small 
region (high dust-to-gas ratio and high frequency) modes are expo- 
nentially damped without oscillating. The other three solutions Fig. 2 
involve linear growth. 

3. THE MODEL 

We work primarily in the thin disk approximation, us- 
ing the vertically integrated equations of hydrodynam- 
ics 
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In these equations, Sg and 17^ are the vertically inte- 
grated gas and dust densities, respectively; u stands for 
the velocity of the gas parcels, P is the pressure, and is 
the gravitational potential. S = Cy ih\P -'yXnSg) is the 
gas entropy, where Cy is the specific heat at constant vol- 
ume and 7 = Cp /Co is the adiabatic index, with Cp the 
heat capacity at constant pressure. T stands for the gas 
temperature. 
The dust evolves Lagrangianly according to 
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where x is the position of a dust particle and v its 
velocity. The gravitational potential is given by ^ = 
-GM*/r^, where G is the gravitational constant. Mi, the 
stellar mass, and r the stellocentric distance. For the 
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Fig. 2. — The three of the five solutions of Eq. (32) that show linear growth. Growth is restricted to the low dust-to-gas ratio (£ < 1), high 
wavenumber {n > 1) region of parameter space. The growing modes in the solutions shown in the middle and lower panels have non-zero 
associated oscillation frequencies - an overstability. Conversely, the solution shown in the upper panel is a true instability. 



pressure we use the ideal gas law P — SgC^ / j, where 
Cg is the sound speed. 

The model is closed by specifying the drag force by 
which gas and dust interact; and Tp, a simple prescrip- 
tion for the gas temperature set by photoelectric heating. 
These are given by 
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The quantities and Xj, are the dynamical and ther- 
mal coupling times between gas and dust, respectively. 
They have a radial variance to match the Keplerian rate 



Tf = Tfo f2o/f2 
T =Tto ^o/f^ 



(8) 
(9) 



where i? = \/ GM* / is the Keplerian angular fre- 
quency. The quantities Tjtq and tjq are free parameters 
of the model. 

Given that the thermal time is sometimes expected to 
be very low (10^^ orbital periods as estimated in Sect. 2), 
we also run models with instantaneous thermal cou- 
pling. For these models, we skip solving the energy 
equation, and equate T = Tp according to Eq. (7). The 



sound speed is updated accordingly. This change effec- 
tively amounts to choosing a new equation of state that 
depends on the dust density. 



lim P = c^{J-l)ToSgSd/So■ 
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We solve the equations with the PENCIL CODE ^ 
which integrates the evolution equations with sixth or- 
der spatial derivatives, and a third order Runge-Kutta 
time integrator. Sixth-order hyperdissipation terms are 
added to Eq. (l)-Eq. (3), to provide extra dissipation 
near the grid scale, explained in Lyra et al. (2008, 2009). 
They are needed because the high order scheme of the 
Pencil Code has little overall numerical dissipation (Mc- 
Nallyet al. 2012). 

4. LINEAR STABILITY ANALYSIS 

We start by performing a linear stability analysis, that 
should assist on interpreting the results of the numer- 
ical simulations. To derive the perturbation equations, 
we make use of the shearing sheet and fluid approxima- 
tions. The first treats the equations in a local, co-rotating 
Cartesian frame. The second greatly simplifies the treat- 
ment of solid particles by having a continuity equation. 
The 2D equations are 

^ The code, including improvements done for the present work, is 
publicly available under a GNU open source license and can be down- 
loaded at http://www.nordita.org/ software/ pencil-code 
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where £ = S^j/S^ is the dust-to-gas ratio and 

Viv = dt + w - V - qQxdy 

is the shear-modified advective derivative, with q = 
3/2 the Keplerian shear rate. Upon linear decomposi- 
tion ip = ipo + ip' and considering axis-symmetric pla- 
nar wave perturbations ip' = i^exp(sf + ikx), these equa- 
tions become 



1 

SV:, = 2nVy (Vx-U^) 



SVy 



1 1 

-iz^Vx (^y-"y) 

2 "^7 

-Erroikux 



(18) 
(19) 

(20) 
(21) 
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where we used the instantaneous thermal coupling ap- 
proximation (Eq. 10) to substitute 



VP = C{EgVEa + EaVEg), 
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with C = c„{y-l)To/Eo = c^f^/ijEo). Equations (18) 
and (21) readily allow for reducing the system to only 
four equations. We substitute these in the radial equa- 
tion for gas velocity to obtain 
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where we also substituted 
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Fig. 3. — Using x — en^ and taking the limit j » 1 allows for better 
visualizing the three behaviours: true instability (red), overstability 
(dark yellow), and damped oscillations (blue). The other two solu- 
tions are the complex conjugate of the oscillating solutions, and not 
shown. 
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We now substitute = 1/^2 so that the dispersion 
relation becomes simpler yet still captures the physi- 
cally interesting case of the most mobile dust. Non- 
dimensionalizing oj = s/ il, n = N / f2 we solve the 
eigenvalue problem 

(M-Ia;)-£r = (30) 
where cr = {vx,Vy,Ux,UyY , I is the unit matrix, and 
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The dispersion relation for this linear system is a quin- 
tic polynomial 

Acv^ + Ba;4 + Cur" + Dup- + Eo; + F = (32) 
with coefficients 



A = X 
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B = 2£ + 2, 




(34) 


C = £^+£{n^ + 2] 
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E = £^{2rp- + \) + 
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r 2 2 2 

F = £ n - £n . 




(38) 



We explore the solutions of the system in the next sec- 
tions. 

4.1. Growth 

Eq. (32) has five solutions. Since qurntics do not have a 
complete analytical solution, we solve it numerically at 
first, to explore the behaviour of the solutions. Since we 
considered planar wave modes of form exp(sf + ikx), 
the imaginary part of the complex root means oscilla- 
tions, whereas the real positive (negative) part translates 
into exponential growth (damping). 
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Fig. 4. — The analytical prediction (Eq. 41) of the linear instabil- 
ity growth compared to the growth rates measured numerically. The 
overall agreement is excellent. The growth rates are only very slightly 
underestimated. 
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Fig. 5. — Although there is no linear instability for e = l, growth oc- 
curs when the amplitude of the initial perturbation is increased. The 
saturation level also depends on the strength of the initial perturba- 
tion. These are hallmarks of nonlinear instability. 

Two of the solutions do not have unstable linear 
modes, representing damped oscillations (Fig. 1). We 
plot in Fig. 2 the three solutions that do show linear 
growth, as a function of the dust-to-gas ratio e and the 
normalized wavenumber n. The left panels show the 
growth rates, the middle ones the damping rates, and 
the right panels the oscillation frequency of the modes. 

The solutions show that there is no linear growth for 
£ > 1. No growth occurs for n < 1 either The linear 
instability is restricted to the phase space n > 1, £ < 1 
only. At low dust load and high wavenumber the three 
growing modes appear. The growing modes in the up- 
per panel have zero oscillation frequency, characteriz- 
ing a true instability As for the two other growing so- 
lutions, the associated non-zero oscillation frequencies 
make them overstabilities. The pattern of larger growth 
rates at large n and low £ invites to explore the limit 
n ^ 1 and £ ^ 1. We thus take x = erp- as characteris- 
tic variable, and explore the behaviour of x ^ 1. In this 



approximation, the coefficients are 

A = l; B=2; C = x; 
D = 3x; E = 3x; F = -x. (39) 

The solutions are plotted in Fig. 3. The instability is 
plotted as the red line. The growth rate is i7/4 for all 
X. The overstability is plotted as the dark yellow line. 
It reaches an asymptotic growth rate of J7/2, at ever 
growing oscillation frequencies. Damped oscillations at 
the epicyclic frequency show as the blue line. The other 
two roots are the complex conjugate of the damped and 
overstable solutions, and not shown. 

4.2. Long wavelength limit - undamped oscillations 

Of particular interest are the free, undamped, oscilla- 
tion that occur through most of the parameter space of 
Fig. 2. We can find these modes analytically by taking 
the long wavelength limit {n ^ 0). Ignoring the n terms, 
the system becomes a quartic equation, 

A = l; B = 2£ + 2; C = £^ + 2£ + 3; 

D = 2£ + 2; E = £^ + 2£ + 2; F = 0; 
which can be solved exactly. The roots are 

lun CO = ±i; lim cc = -(£ -|- 1) ± f (40) 

i.e., two solutions are undamped oscillations at wave 
frequency f2, thus constituting epicyclic oscillations. 
The other two solutions are damped oscillations at the 
same frequency i7, and damping time (e + l)i7~^. The 
eigenvector corresponding to the = ±1 solution rep- 
resent the particular mode for which Vx = Ux and Vy = 
Uy, therefore canceling the drag force. This epicyclic 
mode is low-A: and undamped, and should show up in 
the simulations as a standing wave. We will go back to 
this mode in the next sections. 

5. NONLINEAR NUMERICAL SIMULATIONS 
5.1. Comparing linear theory and simulations 

From the solutions, we see that there is significant 
growth even for very small wavenumbers. The simu- 
lations, however, will cap power at and near the grid 
scale. To make for a meaningful comparison, we add ar- 
tificial Laplacian viscosity v to the gas and dust momen- 
tum equations. The extra term in Fourier space is pro- 
portional to vk'^, which, using the alpha-viscosity recipe 
of Shakura & Sunyaev (1963) v = OiCgH, normalizing by 
Q, and substituting Eq. (29), reduces io q = ajn^. This 
enters in the coefficient matrix as diagonal terms. The 
new, viscous, system is therefore 

[M-l{co + q)]-(r = (41) 

We set a = 10"^ and solve the system numerically. A 
comparison between the predicted and measured linear 
growth rates is shown in Fig. 4. The agreement is excel- 
lent, with the measured growth rates only very slightly 
systematically offset from the analytical prediction. 
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Fig. 6. — End state (left panels) and time evolution (right panels) of the heating instability in one dimension. The photoelectric heating associated 
with the dust generates pressure maxima at the location of dust, that in turn attracts more dust via the drag force. The heating of the gas leads to 
its expansion; the instability maintains a spatial segregation between the gas and the dust. The model shown has t^o = 1, for the most mobile dust, 
and Tjo = 2n. 



5.2. Nonlinear growth 

As seen from Figs. 2 and 4, there is no linear growth 
at e = 1 and beyond, so any growth in this regime must 
be nonlinear. Indeed, nonlinear growth is seen to occur 

for £ = 1 when we raise the noise amplitude from 10"^ 
to 10"^ (Fig- 5). The growth rate and saturation value in- 
crease with increasing amplitude of the initial perturba- 
tion, as expected for subcritical instabilities (e.g. Lesur 
& Papaloizou 2010; Lyra & Klahr 2011). 

5.3. Numerical results 

We model a two-dimensional disk on a uniformly 
spaced grid in cylindrical coordinates {r,(p), ranging 
r=[0.5,1.5]ro and (j)=[-TZ,Tz]. We run models with reso- 
lution [N,.,N^]=[256,256]. The units are such that 



GM* = ro = l2o = ^0 = Cp = 1 



(42) 



To isolate the effect of heating, we shut down migra- 
tion of particles by using initially constant density and 
temperature. The sound speed is set at Cso=0.1, and 
the adiabatic index 7 = 1.4. The reference temperature 



is To = c^q/ [cp(7 - 1)]. Substituting the values, we get 

To = 2.5 X 10"^ in code units. 

We use numerical particles to represent the dust phase 
of the disk, with initially ten particles per grid point, so 
a total of Np = 10 x 256^ = 655360 numerical particles. 
These are initially placed in circular Keplerian orbits. 
The dust-to-gas ratio is set to unity. 

We use reflective boundaries with a buffer zone of 
width 0.1 in each radial border that drives the quantities 
to the initial condition at a dynamical timescale. We use 
outflow no-inlet boundary conditions for the radial ve- 
locity, constant gradient for the azimuthal velocity, and 
zero gradient for density and entropy. The simulations 
are usually run for one hundred orbits at rg. 

5.3.1. The instability in one dimension 

To understand the nature of the instability, it is in- 
structive to consider it first in one-dimensional models, 
since it allows a more through exploration of the param- 
eter space. We also shut down the drag force backreac- 
tion, in order to better isolate the effect of heating (but 
as we show later, this term will become of paramount 
importance in the phenomenology of debris disks in the 
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Fig. 7. — Parameter space exploration for the thermal coupling time Tjg in models with Tjq=1 and without backreaction. The pressure build-up 
associated with too short thermal coupling time leads to modification of the centrifugal balance experienced by a gas parcel. In the case of t^=10 
the pressure buildup leads to supersonic motion and consequent disruption of the dust fingers into a turbulent pattern. The models with larger 
TjQ are better behaved as the buildup is slow. After long enough times (not shown) they too develop shock waves that disturb the dust streams. 
The dynamics is considerably different in models including back-reaction (Fig. 9). 
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Fig. 8. — Parameter space exploration for the friction time t^o in models with Tro=10 and without backreaction. There is an apparent symmetry 
with respect to Tj = 1. As the dust decouples (larger Tj) the particles take longer times to move to the local pressure maxima. As the dust couples 
more strongly (lower Tj), a grain longer times to move away from the density maxima (pressure minima) into the pressure maxima. The velocities 
never get too close to sonic so turbulent disruption as in the case of short thermal coupling time does not occur. 



presence of gas, as shown later). 

In Fig. 6 we show the radial profile of different quanti- 
ties by the end of the simulation at 100 orbits (left panel), 
and their time evolution (right panel). This simulation 
has Tj = 10 and Tjr = 1, i.e., a thermal coupling time 
of 10 orbits and a dynamical coupling time of 1 orbit. 
The initial random overdensities quickly grow into mas- 
sive accumulations of dust. The dust heats up the gas, 
which is seen in the correlation between dust density 
and temperature. As a result of the heating, the gas 
expands, seen in the anti-correlation between gas den- 
sity and temperature. Interestingly, this implies an anti- 
correlation between gas and dust densities. This behav- 



ior has been predicted by Besla & Wu (2007). Rings or 
clumps of photoelectric emitting dust should thus be 
relatively gas-free. Conversely, dense rings and blobs 
of gas should be cold and relatively dust-free. 

This highlights a difference between the present anal- 
ysis and that of Klahr & Lrn (2005). As those authors 
considered a time-independent gas density profile, the 
gas does not expand (compress) upon heating (cooling). 
Let us analyze if any qualitative difference should be ex- 
pected due to expansion or compression. 

5.3.2. Ejfect of pdV work 



8 



Lyra & Kuchner 




Fig. 9. — The effect of the backreaction of the drag force is illustrated 
above. The upper panels show the space-time variation of the dust 
density in models of constant friction time and varying thermal time. 
The lower panels show models of varying friction time and constant 
thermal time. Compared to the same models without backreaction 
(Fig. 7 and Fig. 8), we see that the dust streams are better shepherded. 
This is because when radial velocities are excited, the drag force back- 
reaction damps them back to zero at an e-f olding time Tj Eg/ E^f . The 
shepherding of the dust stream should occur as long as this quantity 
is smaller or similar to . 
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Fig. 10. — Model with instantaneous thermal coupling between gas 
and dust. The gas attains much denser concentrations, and the dust 
executes high frequency oscillations. At later times the dust streams 
in the outer disk break into a roughly uniform sheet of dust extending 
from r=1.4 to r=2.0 

The expansion occurs at the speed of sound, the char- 
acteristic timescale being therefore the sound crossing 
time. Because the structures formed are of size H (where 
H = Cg / f2 is the scale height), the sound crossing time 
equals the dynamic timescale T^y^ = 1/^2. We have 
therefore two timescales; the time scale of the thermal 



coupling Tj, and the timescale of the response of the 
gas, l/i7. Their ratio Tj-ZT^yn = T^i7 is a dimension- 
less number between the second and first terms in the 
right hand side of the entropy equation (Eq. 3), namely 
relaxation and advection. Hidden in the latter term is 
the pdV work. If f2 is too big, relaxation is slow and 
advection dominates, i.e., the expansion is adiabatic. If 
on the other hand t^-]? is too small, advection is slow 
and relaxation dominates. That means that the temper- 
ature changes faster than the gas can respond. The gas 
is unchanged and the process is therefore isochoric. So, 
schematically, the temperature rise from T2 to Ti will, in 
each regime, brings the pressure from Pi to P2 according 
to 

if T^/7»l =^ P2 = J'i(T2/Ti)')'/(t-i) (43) 
if T^/2<1 ^ P2 = Pi{T2/Ti) (44) 

and as 7 > 1, P2/P1 arid T2/T1 have the same sign, so a 
rise in temperature leads to a rise in pressure irrespec- 
tively of the density drop. The instability would only be 
quenched if the temperature rise occurred at constant 
pressure. This condition could only be fulfilled if the 
expansion occurred in a huge pressure reservoir, which 
is not the case of gas-poor debris disks. We are safe to 
state that the photoelectric heating of the gas will lead 
to pressure rise and further concentration of dust if the 
latter is mobile. We show this in Fig. 7, where models 
with the exploration of the parameter space for T^f2 is 
shown. From left to right we show models of t^- i7=l, 10, 
1000, and 1000. The friction time was kept constant at 
Ty-i7=l so that the dust is mobile. In all cases the dust 
concentrates. 

The models shown in Fig. 7 contain some interest- 
ing features worth highlighting. As seen in the model 
of T^i7=1.0, there are instances in time, around 30 or- 
bits, that the dust distribution rapidly migrates inwards, 
setting in another equilibrium location. In the model 
with Tj]7=10 the dust at later times, around 50 orbits, 
passes from a thin stream to a shower, resembling the 
transition to turbulence seeing in, e.g., cigarette smoke. 
Both effects seem to be related to the pressure buildup. 
As the dust concentrates, the local temperature rises, 
leading to further dust buildup. The runaway process 
has to saturate at some point, and these effects seem to 
be manifestations of saturation. The pressure buildup 
leads to a gradual change in the centrifugal equilib- 
rium 

= rill + dyP. At some point the buildup 
of pressure significantly changes that relation, and the 
disk readjusts. In the other case, the buildup of pres- 
sure leads to supersonic speeds. Shock waves develop, 
leading to the disruption of the dust stream. These ef- 
fects will be significantly mitigated by the inclusion of 
the drag force backreaction, in Sect. 5.3.4. 

5.3.3. Effect of gas drag 

The models shown in the previous section were run 
with a fixed friction time of TfQ=l, representing the 
most mobile particles. We consider now models with 
varying Tj, keeping Tj-i7=10. A suite of such models is 
shown in Fig. 8. As the dust decouples, it takes time 
for the pressure gradient to move the particles to the 
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Fig. 11. — Snapshots of the dust density, at 10, 50, and 200 orbits, respectively. The heating-clumping instability concentrates the dust axis- 
symmetrically into rings. The axis-symmetric is due to strong damping of non-axissymmetric modes provided by the drag force backreaction. 
At latter times, one of these rings is seen to oscillate, appearing as an eccentric ring. This is not because of the presence of a eccentric planet 
shepherding the ring, but simply the result of a standing wave propagating along it. 



pressure maxima. As the dust couples more strongly, 
it takes time for the dust to decouple itself from the den- 
sity maxima (pressure minima) and into the pressure 
maxima. The symmetry with respect to T^i7 = 1 is strik- 
ing. 

In the extreme of particle tracers (t^ = 0) there is no 
instability as the gas and dust cannot get separated. The 
heating leads to expansion, that carries away the dust, 
leading to cooling. The process is self -regulated. In the 
opposite extreme of decoupled dust ^ oo there is no 
instability either as the dust does not get pushed toward 
pressure maxima. 

5.3.4. Ejfect of gas drag back-reaction 

We now turn to the effect of the last term in the right 
hand side of the momentum equation, the drag force 
backreaction from the particles onto the gas. We re-run 
models of T-j^f2=l and 10, with TjrQ=l (the leftmost ones 
in Fig. 7), but now including this term. The results are 
shown in the upper panels of Fig. 9, for the dust den- 
sity only. Comparing these models with the ones with- 
out back-reaction shown in Fig. 7, we see that the jerks 
in position are absent, as well as the dispersion of the 
particles into turbulent streams as seen in the models 
without backreaction. In other words, the backreaction 
of the drag force has the effect of shepherding the dust 
streams. 

The reason is because in the other models, although 
the dust is forced to follow the gas, the gas is uncon- 
strained by the dust. When backreaction is included, 
the gas follows the dust whenever the dust-to-gas ra- 
tio £ is high. When £ ^ 1 the backreaction dominating 
the gas motion effectively herds the dust; a linear per- 
turbation to the system executes exponentially damped 
oscillations, as discussed in Sect. 4. 

Interestingly, this effect, being proportional to t^, 
should break the symmetry with respect to = 1 seen 
in the models without backreaction (Fig. 8). Indeed, this 
is what is seen in the lower panels of Fig. 9, where we 
show the result of two simulations with t^=10 and back- 
reaction included, varying t^. The left panel has tj = 0.1 



and the right panel has Tjr = 10. The simulation with 
shorter friction time experiences less clumping for the 
same thermal coupling time. 

The herding provided by the drag force also intro- 
duces a dependency on Tj. / Tf, because if <^ Tf, the 
pressure builds much faster than the dust can respond. 
The backreaction becomes less important, and the situa- 
tion becomes similar to the one shown in the left panels 
of Fig. 7. 

5.3.5. A model with instantaneous thermal coupling 

We test now models with instantaneous thermal cou- 
pling, where the equation of state is given by Eq. (10). 
The time evolution of the dust and gas densities is 
shown in Fig. 10. The evolution looks a hybrid be- 
tween the models shown so far On the one hand the 
dust streams are very thin, because of the shepherd- 
ing provided by the drag force backreaction. On the 
other hand, the infinitesimally small thermal time leads 
to rapid wave excitation, seen in the high frequency os- 
cillation of the dust. Interestingly, at later times, (<80 
orbits) the thin streams in the outer disk get disrupted 
and form a roughly uniform extended sheet of particles 
spanning r = [1.4,2.0], a significant part of the radial do- 
main. A very dense concentration of gas {Sg = 60) oc- 
curs in the low-pressure dust void left by the massive 
particle concentration at r=1.2 

5.3.6. The model in two dimensions 

Having analyzed the behavior of the instability in ID 
models, we turn our attention to two-dimensional mod- 
els. We show in Fig. 11 snapshots of the disk at 10, 50, 
and 185 orbits. The appearance of the disk in Carte- 
sian projection, for the last snapshot, is shown in Fig. 12. 
Though the initial distribution is clumpy, the system re- 
arranges itself in rings, which suggests a damping of 
non-axissymmetric modes. We checked that in the ab- 
sence of backreaction there is no development of rings 
- the system consists simply of clumps. This implies 
that it is the drag force backreaction that is providing 
the damping. The regular spacing of the rings, also seen 
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Fig. 12. — Appearance of the disk in Cartesian coordinates. An oscillation traveling along the ring between r=0. 8-0.9 shows itself as deviations 
from axis-symmetry, at 4 and 11 o'clock. The eccentricity of the ring is e = 0.04, similar to the eccentricity observed in the debris disk around 
HD 61005. 



in the one-dimensional models, speak of a preferred ra- 
dial wavelength for the instability (though no hint of 
this most unstable wavelength was found in the linear 
stability analysis). 

The most interesting development of the model, 
though, is seen in the rightmost panel of Fig. 11. At 
late times, one of the rings starts to oscillate. The os- 
cillation does not grow in time, which shows that it is 
not an instability. In fact, as shown before, it is rather 
an epicycle executed in the orbital plane, with a period 
equaling the Keplerian period. The observational up- 
shot is that the disk appears off-centered. In Fig. 13 we 
show this particular ring only. The location of the star 
[x,y] — [0,0] is marked with a red "x". The dusty ansae 
appear offcentered, the lower one at minimum elonga- 
tion rinin=0.84, and the upper one at maximum elonga- 
tion rmax=0.91. These "apsides" are not exactly diamet- 
rically opposed. Though they are extended, we define 



their azimuthal location as the azimuths of the points 
of minimum and maximum elongation, (pmax and (^min/ 
respectively. The "center" {rc,(pc\ of the eccentric ring 
is thus at rc=rmdLX - J'min/2, the middle point in radius; 
and (pc = (pmax + (tt- ^min)/2, the average between the 
azimuth of the outer ansa and the azimuth of (the an- 
tipode) of the inner ansa. This is marked as a black "x". 
The eccentricity of the ring is 



= 0.04. 



(45) 



a value that matches the eccentricity found for the ring 
around HD 61005 (e=0.045 ± 0.015, Buenzli et al. 2010. 
Though they fall short of the eccentricity of the Foma- 
Ihaut ring (e = 0.11, Kalas et al. 2005), larger ampli- 
tude oscillations could account for that eccentricity. The 
point to make is that a standing wave through a narrow 
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Fig. 13. — Appearance of the oscillating ring. The star and the center 
of the ring defined by the dusty ansae are marked red and black "x", 
respectively. Due to the wave propagating along it, the ring appears 
off-centered, with eccentricity e = 0.04. 



Dust density 




Indeed, this situation could explain the bright moving 
source in the debris disk around Fomalhaut. The source 
was initially thought to be a planet shepherding an ec- 
centric sharp ring (Kalas et al. 2005), but that hypothesis 
was dealt a serious blow recently as the source could 
not be detected in the near infrared 0anson et al. 2012). 
The most likely explanation, as suggested by the au- 
thors, is an extended cold dust cloud shining through 
reflected light. Recent work by (Boley et al. 2012) sug- 
gests that the ring is confined by a pair of shepherding 
terrestrial mass planets, well below the current detec- 
tion limits. Detection of gas around the ring would be 
a way to distinguish that scenario from the one we pro- 
pose. At present, only upper limits on the amount of gas 
in the Fomalhaut system exist (Liseau 1999); they are rel- 
atively insensitive, however, because they probed CO 
emission, and CO could easily be dissociated around 
this early A type star. 

We notice also that at times the ring displays multi- 
ple apsides, which means oscillations with higher fre- 
quency than simple epicycles. Undamped modes with 
ci? « A: appeared in the linear analysis for £ = 1 and, 
though the dust-to-gas ratio is not strictly unity in this 
d5mamical case, similar undamped modes should ac- 
count for this behavior. 

We modeled also a r - z slice of the disk, in order to 
investigate the vertical structure of the debris disk. The 
resolution is 256x128, spanning 2H above and below 
the midplane, and 8H in the radial direction. Particles 
of three sizes were used, of friction times T|=0.0,1 0.1 
and 1.0. The behavior of the particles resembles that in 
the simulations of Johansen et al. (2006) for particle set- 
tling in the midplane of protoplanetary disks. The more 
decoupled particles settle first, and stir the midplane 
into a mild Kelvin-Helmholtz turbulence that supports 
through diffusion a particle subdisk of finite thickness 
(Dubrulle et al. 1995; Garaud & Lin 2004). The differ- 
ence in the debris disk case reflects the fact that the dust 
is the heating agent, so once the dust settles, the gas in 
the midplane starts to get heated up. This gas expands 
and eventually becomes buoyant. At later times, as seen 
in the rightmost panel of Fig. 15, the midplane of the 
disk develops significant voids of gas, that is pressed 
against the upper layers. This may account for the dis- 
parity found in the measured gas scale heights in the 
few debris disks where gas lines could be measured. 



Fig. 14. — The maximum dust density in the model shown in Fig. 11. 
The peak near the end of the simulation is the bright dust clump lo- 
cated at »i 1.2 seen in the rightmost panel of Fig. 11. In actual images 
of debris disks, this clump could be bright enough to be mistaken for 
a planet shepherding a sharp eccentric ring. 



ring is seen as eccentricity, a feature which is usually 
attributed to the gravitational influence of hjrpothetical 
undetected planets in eccentric orbits. 

We also see that although axis-symmetry is strongly 
favored, some rings broke into clumps at later times. 
In the rightmost panel, at t = 185 orbits, the clumps at 
r ~ 1.2 achieved a very high dust load with S^/Sg ~ 40, 
sometimes exceeding 100 (see Fig. 14). This is interest- 
ing because these dust clumps are seen quite close to 
the borders of the narrow rings, and could be mistaken 
for a planet due to their brightness and compactness. 



6. CONCLUSIONS 

We modeled for the first time the hydrodynamics of 
debris disks. These disks are usually assumed to be gas- 
free, yet outgassing from planetesimals and dust grains 
should occur, supplying gas and maintaining a dust-to- 
gas ratio of roughly unity. 

Heated by the star, the dust transfers the heat to 
the gas, which expands, contracting again in colder re- 
gions devoid of dust. An anti-correlation between dust 
and gas densities follows, which was not observed by 
(Klahr & Lin 2005). In one dimension, the system de- 
velops well-spaced spikes in dust density, that in two- 
dimensions assume the shape of well-defined, sharp 
rings. The sharpness is a consequence of the drag force 
backreaction, that efficiently damps high frequency os- 
cillations at high dust load. This herding by the drag 
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Fig. 15. — A simulation in the r-z plane, in order to investigate the vertical structure of debris disks subject to the heating-clumping instability. 
The upper panels show the gas density, and the lower the dust density. Particles of three sizes are represented. The dust settles to the mid- 
plane, supported by self-sustained Kelvin-Helmholtz turbulence. It heats the gas in the midplane, that in turn expands. At later times, dust-gas 
segregation has been established, with the hot dust expelling the gas from the midplane. 




Fig. 16. — A snapshot of the model, inclined to the same orientation 
of the debris disks around Fomalhaut. The bright dots are high con- 
centrations of dust, and would be very bright in reflected light. As 
they also move in Keplerian orbits, they could be easily mistaken for 
planets. Scan the barcode at the top left corner for an animation of the 
simulation. 

force has not previously been accounted for, and may 
explain the sharp rings seen in actual images of debris 
disks. A snapshot of the model, rotated to the same in- 
clination of the debris disk aroimd Fomalhaut, is shown 
in Fig. 16. 

The system supports a spectrum of oscillations, most 
of them damped by the drag force. However, some spe- 
cific modes for which the gas and dust velocity coin- 
cide correspond to undamped oscillations. For these 
modes, the rings behave as wave guides. The visual re- 
sult of these oscillations is a finite eccentricity. In our 
model, we measure e = 0.04 for one of the oscillatrng 
rings. This is in good agreements with the observed 
eccentricities of rings in debris disks. These eccentric- 
ities are usually presented as convincing evidence for 



the presence of unseen planets in eccentric orbits, whose 
gravitational influence would shepherd the ring and im- 
part its eccentricity to it. Our model explains both fea- 
tures without the need of a planet. The eccentricity 
is simply a standing wave propagating along the ring. 
This fits neatly with the recent (Janson et al. 2012) in- 
frared non-detection from the bright source around Fo- 
malhaut (Kalas et al. 2005), indicating that the source 
is most likely a dust cloud and not a planet shepherd- 
ing the eccentric ring, as previously assumed. We notice 
that the possibility of these features arising from density 
waves instead of from gravitational influence of unseen 
companions has been raised by Jalali & Tremaine (2012). 
The difference between that work and the present one is 
that they ignored dust-gas interactions, relying on stel- 
lar encounters to trigger the density waves. We show 
instead that oscillations can arise through hydrodynam- 
ics, without the need for external perturbers. 

We also model a r - z slice of the disk, following the 
settling of the dust and how it impacts the gas. We find 
that the gas is expelled from the midplane. It presses 
against the gas above and below, increasing the gas 
density in the atmosphere. The vertical structure of 
debris disks is therefore hot thin gas in the midplane, 
cold dense gas above and below, as also previously sug- 
gested by Besla & Wu (2007). This highly unstable situ- 
ation is maintained as long as the dust is present in the 
midplane. 

We caution that the models presented are very sim- 
plified, and should be taken more as proofs of concepts 
than conclusive evidence. We plan in future publica- 
tions to include magnetic fields, add a more realistic 
treatment of the disk thermodynamics, perform a full 
3D analytical analysis and numerical calculations, and 
push the resolution to the limit of modern computa- 
tional power. Still, that the simple models presented 
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here shed new light in pertinent unsolved questions of 
debris disks is an illustration of the importance of hy- 
drodynamics in these systems. 
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